/*
 * EquilibriumProblem.h
 *
 *  Created on: 27 Jul 2011
 *      Author: Allan
 */

#ifndef EQUILIBRIUMPROBLEM_H_
#define EQUILIBRIUMPROBLEM_H_

// Eigen includes
#include <Eigen/Core>
using namespace Eigen;

// GeoReact includes
#include <Numerics/NewtonSolver.h>

class EquilibriumProblem : public NewtonProblem
{
public:
	EquilibriumProblem(const MatrixXd& Aii,
		const MatrixXd& uei,
		const MatrixXd& alpha,
		const VectorXd& hj,
		const VectorXd& ln_a,
		const VectorXd& ln_Ke);
	
	virtual	~EquilibriumProblem();
	
	virtual void
	Function(const VectorXd& n, VectorXd& F);
	
	virtual void
	Jacobian(const VectorXd& n, MatrixXd& J);
	
private:
	/// The number of species (N), equilibrium reactions (Ne) and equilibrium conditions (Nj)
	unsigned N, Ne, Nj;
	
	/// The matrix of activity derivatives of the species
	const MatrixXd& Aii;
	
	/// The stoichiometric matrix of the system of reactions
	const MatrixXd& uei;
	
	/// The matrices that define the set of equilibrium constraints in the system 
	const MatrixXd& alpha;
	
	/// The right-hand side vector hj of the system of equilibrium conditions
	const VectorXd& hj;
	
	/// The natural log of the activities of all species
	const VectorXd& ln_a;
	
	/// The natural log of the equilibrium constants of the equilibrium-controlled reactions
	const VectorXd& ln_Ke;
};

#endif /* EQUILIBRIUMPROBLEM_H_ */
